capture log close

clear all

cd ${path}analysis
set more off
set matsize 2000

log using smog_pollution_pm10.log, replace


use pollution_star_merged

xtset county date

/*First, divide the counts of inspections by 1000, to make the coefficients easier to read*/

gen lag = 0
foreach var of varlist firstgen*  othergen* {
	qui replace `var' = (`var')/1000
}


/* Generate day-of-week and calendar week*/

gen dow = dow(date)
gen week = wofd(date)



/*Declare the "if" conditions in a local, as well as the full set of controls*/

local conditions  if year>=1998 & inlist(county,5,12,18,22,26,53)==0 & (inlist(county,3,6,11,12,13,17,23,47,58)==0 | (year !=2002&year!=2003)) 
local controls tMin tMax prec 
local absorbvars i.county i.week i.county#i.year 

/*PM10 is measured weekly-ish.  Take the average reading and weather, and the 
count of re-inspections from the last PM10 reading (some counties have more than 1 per 
week because they have more than 1 sensor)*/

keep if pm10 !=.

gcollapse pm10 tMin tMax prec (last) firstgen* othergen* fpr*, by(county year week)


/*Columns to store bin coefficients*/
gen newcars = _n>20 in 1/40
gen bin = _n-1 -20*(_n>20) in 1/40

/*Base regression of re-inspections*/

	eststo pm10_base: reghdfe pm10 firstgen_90 othergen_90 `controls' `conditions', vce(cluster county) absorb(`absorbvars')
	
	test firstgen_90 = othergen_90
	estadd scalar oldnew_F = r(F) : pm10_base
	estadd scalar oldnew_p = r(p) : pm10_base
	
/*Interaction with station_quality*/

eststo pm10_fpr_interaction_base: reghdfe pm10 c.fpr_current_firstgen_90#c.firstgen_90 firstgen_90 c.fpr_current_othergen_90#c.othergen_90 othergen_90 `controls'  `conditions', vce(cluster county) absorb(`absorbvars')
	/* Calculate marginal effect at the mean */
	margins, dydx(firstgen othergen) at((mean) fpr_current_firstgen_90) noestimcheck  post

	foreach gen in firstgen othergen{
		local `gen'_ame = _b[`gen'_90]
		local `gen'_se=_se[`gen'_90]
		estadd scalar `gen'_ame = ``gen'_ame' : pm10_fpr_interaction_base
		estadd scalar `gen'_se = ``gen'_se' : pm10_fpr_interaction_base
		global pm10_`gen'_star = cond(abs(_b[`gen'_90]/_se[`gen'_90]) < invt(e(N),.95),"",cond(abs(_b[`gen'_90]/_se[`gen'_90]) < invt(e(N),.975),"*",cond(abs(_b[`gen'_90]/_se[`gen'_90]) < invt(e(N),.995),"**","***")))
	
	}

	test firstgen_90 = othergen_90
	estadd scalar oldnew_F = r(chi2) : pm10_fpr_interaction_base
	estadd scalar oldnew_p = r(p) : pm10_fpr_interaction_base	
/*Quality bins*/

	eststo pm10_bin: regress pm10 firstgen_fprbin*_90 othergen_fprbin*_90 `controls' `absorbvars'  `conditions', vce(cluster county)
	
	gen pm10_fprbin_coef=.
	gen pm10_fprbin_ub = .
	gen pm10_fprbin_lb = .
	local row = 1
	forvalues i = 0/19{
		replace pm10_fprbin_coef = _b[firstgen_fprbin`i'_90] if newcar == 0 & bin == `i'
		replace pm10_fprbin_coef = _b[othergen_fprbin`i'_90] if newcar == 1 & bin == `i'
		
		replace pm10_fprbin_ub =_b[firstgen_fprbin`i'_90]+ _se[firstgen_fprbin`i'_90] * invnormal(1-.025/20) if newcar == 0 & bin == `i'
		replace pm10_fprbin_lb =_b[firstgen_fprbin`i'_90]- _se[firstgen_fprbin`i'_90] * invnormal(1-.025/20) if newcar == 0 & bin == `i'
		
		replace pm10_fprbin_ub = _b[othergen_fprbin`i'_90] + _se[othergen_fprbin`i'_90] * invnormal(1-.025/20) if newcar == 1 & bin == `i'
		replace pm10_fprbin_lb = _b[othergen_fprbin`i'_90] - _se[othergen_fprbin`i'_90] * invnormal(1-.025/20) if newcar == 1 & bin == `i'

	}
	
	
/* Table 5: PM10 results*/

#delimit ;

esttab pm10_base  pm10_fpr_interaction_base  using pm10.tex, replace se booktabs label nomtitles  star( * .1 ** .05 *** .01)
drop(tMin tMax prec *.county* *.week *date _cons *othergen*, relax ) 
stats(firstgen_ame firstgen_se, label("\emph{Effect at Mean Station Quality}" " ") layout(\emph{@}\sym{${pm10_firstgen_star}} (\emph{@})))
title(Station Quality and County-Level PM10 Pollution\label{tab:pm10})
varlabels(firstgen_90 "000s of Re-Inspections Last 90 Days\\ \quad\quad 1975-1985 Vehicles"  othergen_90 "\quad\quad 1985+ Vehicles"
c.fpr_current_firstgen_90#c.firstgen_90 "\quad\quad 1975-1985 Vehicles \$\cdot\$ Station Quality" 
c.fpr_current_othergen_90#c.othergen_90 "\quad\quad 1985+ Vehicles \$\cdot\$ Station Quality")
order(firstgen_90 c.fpr_current_firstgen_90#c.firstgen_90 ) 
postfoot( ) prefoot( ) nocons noobs mgroup("Outcome is PM\$_{10}\$ (\$\mu/m^3\$)", pattern(1 0 0 0 0) span prefix(\multicolumn{@span}{c}{) suffix(}))
note(Note: Observations are county-weeks.  All regressions control for daily weather, county fixed effects, calendar week fixed effects, and county-year fixed effects. Standard errors clustered by county reported in parentheses.)
substitute( "\begin{tabular}" "\setlength{\linewidth}{.1cm}\newcommand{\contents}{\begin{tabular}"
            "\end{tabular}" "\end{tabular}}\setbox0=\hbox{\contents}\setlength{\linewidth}{\wd0-2\tabcolsep-.25em}\contents"
            "{l}{\footnotesize" "{p{\linewidth}}{\footnotesize")

;

esttab pm10_base  pm10_fpr_interaction_base  using pm10.tex, append se booktabs label nomtitles  star( * .1 ** .05 *** .01)
drop(tMin tMax prec *.county* *.week *date _cons *firstgen*, relax )  
stats(othergen_ame othergen_se oldnew_F oldnew_p, label("\emph{Effect at Mean Station Quality}" " " "F-test of Equality" "\$P > F \$") layout(\emph{@}\sym{${pm10_othergen_star}}  (\emph{@}) @ @))
varlabels(firstgen_90 "000s of Re-Inspections Last 90 Days\\ \quad\quad 1975-1985 Vehicles"  othergen_90 "\quad\quad 1985+ Vehicles"
c.fpr_current_firstgen_90#c.firstgen_90 "\quad\quad 1975-1985 Vehicles \$\cdot\$ Station Quality" 
c.fpr_current_othergen_90#c.othergen_90 "\quad\quad 1985+ Vehicles \$\cdot\$ Station Quality")
order(  othergen_90 c.fpr_current_othergen_90#c.othergen_90)
prehead( ) posthead( ) prefoot( ) nonum  nocons postfoot(\bottomrule @starlegend \\ \multicolumn{@span}{l}{\footnotesize @note} \end{tabular}\end{table})
note(Note: Observations are county-weeks.  All regressions control for daily weather, county fixed effects, calendar week fixed effects, and county-year fixed effects. Standard errors clustered by county reported in parentheses. F-tests of equality report the F-statistic from testing the null hypothesis that either the coefficient (column 1) or the average marginal effect (column 2) of 1975--1985 model year vehicles equals that of 1985+ model year vehicles.)
substitute( "\begin{tabular}" "\setlength{\linewidth}{.1cm}\newcommand{\contents}{\begin{tabular}"
            "\end{tabular}" "\end{tabular}}\setbox0=\hbox{\contents}\setlength{\linewidth}{\wd0-2\tabcolsep-.25em}\contents"
            "{l}{\footnotesize" "{p{\linewidth}}{\footnotesize")
;

replace bin = bin*.05;

/* Figure A7: Bin plot*/

#delimit ;
local mu = char(181);

graph twoway scatter pm10_fprbin_coef bin if newcar == 0 , msize(small) || scatter pm10_fprbin_coef bin if newcar == 1, msize(small) msymbol(S)||
lowess pm10_fprbin_coef bin if newcar == 0 || lowess pm10_fprbin_coef bin if newcar == 1, lpattern(dash)||,
legend(order(- "1975-1985 Vehicles" - "1985+ Vehicles" 1 2 3 4) label(1 "Coefficients") label(2 "Coefficients") label(3 "Lowess Best Fit") label(4 "Lowess Best Fit")) 
ytitle("Change in PM10 ({&mu}/m{superscript:3}) per 1000 Re-inspections")
xtitle(Bottom of Station Quality Bin) graphregion(color(white)) title(PM10) name(nox, replace)
;

graph export pm10_bin.pdf, replace;

